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INTRODUCTION 


A  great  deal  of  Interest  has  been  shown  within  the  Army  in  recent  years  in 
gun  fired  projectiles  which  correct  their  flight  paths  by  impulsive  means.  One 
of  the  more  interesting  concepts  for  large  caliber  projectiles  has  been  an 
arrangement  of  propellant-filled  helical  grooves  on  the  cylindrical  portion  of 
the  projectile  which  burn  in  a  progressive  manner  (from  one  end  to  the  other). 
This  arrangement  has  the  virtue  of  controlling  the  introduction  of  nutation  and 
precession  motion  to  the  flight. 

In  the  flight  simulation  of  such  projectiles,  the  burning  of  the  propellant 
has  historically  been  modeled  as  the  application  of  a  radial  force  to  the  projec¬ 
tile  which  moves  from  one  end  of  the  groove  to  the  other  as  the  point  of  burning 
of  the  propellant  moves.  Since  the  typical  speed  of  burning  might  be  20,000  feet 
per  second,  the  integration  step  size  (in  time)  required  to  adequately  model  the 
shape  of  the  helic  groove  in  space  will  be  quite  short  compared  to  the  typical 
rigid  body  aerodynamic  response  times  of  the  projectile,  causing  very  large  costs 
in  computer  time  in  the  flight  simulation  of  the  projectile.  Therefore,  an 
impulsive  analytical  model  of  such  a  burning  process  was  developed  which  exactly 
models  the  dynamic  response  of  a  rigid  axisymmetric  body  to  which  is  applied  a 
radial  thrust,  moving  from  one  end  to  the  other  with  a  finite  burn  time.  This 
process  is  applied  as  a  single  step  change  in  the  spin  components  in  body  fixed 
coordinates,  resulting  in  large  savings  in  computer  time  compared  to  the  numeri¬ 
cal  integration  of  the  rigid  body  response.  This  is  not  an  approximation  to  the 
rigid  body  motion  und  r  this  force  alone,  but  an  exact  solution;  the  numerical 
solution  is  an  approximation  to  the  geometry  of  the  helix.  If  there  is  an 
approximation,  it  is  in  the  neglect  of  the  aerodynamic  forces  during  the  burn 
time.  Aerodynamic  forces  typically  provide  approximately  two  g's  of  acceleration 
which  is  almost  entirely  axial,  while  the  impulsive  thruster  typically  provides 
three  thousand  radial  g's  during  the  short  duration  of  the  burn.  This  would 
appear  to  amply  justify  neglect  of  the  lateral  aerodynamics  during  the  short 
application  time  of  the  impulse. 


ANALYSIS 


The  Euler  equations  of  rotational  motion  of  a  rigid  body  can  be  written  as 

I  w  -w  w  (I  -I  )  =  N  =  0  (1) 

x  x  y  z  y  z  x 


1  w  -w  w  (I  -I  )  =  N  =  (x  +  Ax)  F  cos(9)  (2) 

y  y  z  x  z  x  y 


Iw-ww(I-I'=N  =  ( x  +  Ax)  F  sin( 9)  (3) 

z  z  x  y  x  y  z 

where  Nx  =  0,  since  the  applied  force  is  radial  and  exerts  no  torque  about  the  x- 
axis,  and  I  -I  =0  since  the  body  is  assumed  axially  symmetric  (fig.  1).  The 
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moment  arm  due_to  the  instantaneous  application  of  the  force  of  the  burning  pro¬ 
pellant  is  at  x  +  Ax,  where  Ax  is  the  displacement  between  the  initiation  point 
of  ignition  and  the  center  of  mass  of  projectile,  and  x  =  0  at  the  point  of  ini¬ 
tiation  of  the  burning  of  the  propellant. 

Several  auxiliary  relationships  are  required  for  this  development  (fig.  2). 
Inspection  of  this  figure  reveals  that 

I  »  t  -  t  (4) 


x  =  x-x,.  +  Ax  =  x-  x 


9  »  (9  -  9Q)  =  2  i  Tx/D 


l2  =  {  x2  +  (rQ)2^ 


_  _  9  1/ 

1  =  -X  [1  +  (irT)  ]'2 


where  1  is_the  distance  down  the  groove,  xQ  is  the  value  of  _x  at  die  start  of  the 

burn,  and  x  is  zero  at  the  start  of  the  burn.  Similarly,  t  and  0  are  both  zero 

at  the  start  of  burn.  The  quantity  xM  is  the  value  for  x  at  the  center  of  mass 

and  Ax  is  the  displacement  in  x  between_the  center  of  mass  and  the  start  of 

thruster  burn.  Note  from  equation  5  that  x  +  Ax  =  x  -  x,.. 

M 

If  the  linear  burning  rate  of  the  propellant  is  R  feet  per  second,  and  it 
yields  I  pound-seconds  of  thrust  per  foot  of  thruster,  then 

F  =  RI  (9) 

is  the  force  of  thrust  being  applied  to  the  groove  at  the  point  of  burning  at  any 
time  during  the  burning  process. 

From  equation  6 


dx/dt  =  —  [ D/ (2i»T)  ]  d9/dt 


From  equation  8 


dl/dt  =  -dx/dt  /  ( 1  +  ( nT)  ) 


n/MVV. 
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dl/dt  =  K 


where  R  is  Che  linear  burning  rate. 


From  equation  1 


dwx  =  0 


for  an  axially  symmetric  projectile, 


From  equation  8 


x  =  -1/  Ml  +  ( rrT)2  ]  =  -Rt/  /  [1  +  (itT)2] 


From  equations  6  and  14 


8  =  9  -  eo  =  +2nTx/D  =  -  2  tt  RT  t  /  {D  /  [1  +  (ttT)  ]} 


Thus,  when  equations  9,  14,  and  15  are  applied,  equations  2  and  3  become 


I  w  -  w  w  (I  -  I  )  =  {-Rt/  /  (i  +  (ttT)  ]  +  Ax}*RI* 
y  y  z  x  z  x 

cos{-2ttTrT/[D*  /(I  +  (ttT)2)]  +  0O} 


I  w  -  w  w  (I  -  I  )  =  {-Rt/  /  {1  +  (wT)  ]  +  Ax}*RI* 
z  z  x  y  x  y 

sin  {-2wTRt/  [  D*  /(I  +  (  xT)  2  )  ]  +  0  q  } 


Define 


A  =  -R2!/  /  (l  +  ( xT) 2  ] 


B  =  R  IA  x 


K  =  -2xTR/  {D*  /  ( 1  +  ( ttT)  )  ] 


then  equations  16  and  17  become 


Iw-ww(I-I)=  [At  +  B]  cos ( Kt  +  9_) 
yyzxzx  0 

Iw-ww(I-I)=  [At  +  B]  sin(Kt  +  9.) 
z  z  y  x  x  y  0 


Equations  21  and  can  be  decoupled  by  differentiating  each  and 
tuting  the  solution  for  w  from  22  into  the  derivative  of  21  and  the 

*  Z  • 

for  w  from  21  into  the  derivative  of  22.  Keeping  in  mind  that  w  = 
differentiation  results  in  x 

Iw-ww(I-I)=A  cos(Kt  +  9,-)  -  AKt  sin(Kt  +  9n) 
yyzxzx  0  0 

-  BK  sin(Kt  +  00) 

Iw-ww(I-I)=A  sin<Kt  +  9_)  +  AKt  cos(Kt  +  0_) 
z  z  y  x  x  y  0  0 

+  BK  cos(Kt  +  eQ) 


The  substitution  of  the  solution  of  22  into  23  eliminates  w 


Iw  +  w  w  (I-I)(I-I)/(I)= 
yy  xyxz  x  y  z 

{At[w (I  -  I)  /  (I  )  -  K)  +  B[w  (I  -  I)  /  (I  )  -  K] } 

A  4*  A  4  A  A  A 


sin(Kt  +  0Q)  +  A  cos(Kt  +  9^) 


while  the  substitution  of  the  solution  of  21  into  24  eliminates  w 


I  w  +  v  w  (I-I)(I-I)/(I)  = 
zzxzxz  xy  y 


(At [ w  (I  -  I  )  /  (I  )  +  K]  +  B[w  (I  -  I  )  /  (I  )  +  K] > 
xxy  y  xxy  y 

cos(Kt  +  Oq)  +  A  sin(Kt  +  9^) 


(21) 

(22) 

substi- 
solution 
0  ,  the 

(23) 

(24) 


(23) 


(26) 


For  convenience,  rewrite  25  as 
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m  =  w  ( r  -  i  )  (i  -  i  )  /  (i  ) 
x  x  z  x  y  z 


E  =  -A[K  +  w  (I  -  I  )  /  (I  )) 
X  X  z  z 


F  =  -B[K  +  wx(Ix-  Iz  )  /  (Iz  )] 


The  complete  solution  to  27  is  the  sum  of  the  general  solution  of  the  equation 
with  the  driving  terms  set  to  zero  and  a  particular  solution  of  the  inhomogeneous 
equation.  The  general  solution  of  the  homogeneous  equation 


Jw  +  Mw  =0 

y  y 


w  =  c,  cos(aC  +  4>) 

y  h 


which,  if  substituted  back  into  the  homogeneous  equation,  yields  the  secular 
equation 


-x  j  +  m  =  n 


X  =  /  (M/J) 


A  particular  solution  of  the  inhomegeneous  equation  is  needed, 
of  the  form 


Try  an  expression 


w  =  cj  sin(Kt  +  0^)  +  c^cos(Kt  +  9^)  +  c^t  sin(Kt  +  9^) 

+  c.T  cos(kT  +  0r ) 

4  0 

Differentiate  35  twice  .o  obtain  w  ,  then  substitute  it  and  35  into  27  to  evalu- 

v 

ate  Cj,  c 2>  c^,  and  c^.  We  get 


•c,K  J  -  2c. KJ  +  Me,  =  F 
1  4  1 


v*-’ 


•HI 


'  *1 


a 


‘.*1 


*1 


>yj 


-c^K  J  +  Mc3  =  E 


-c^K  J  +  2c3KJ  +  c2M  =  A 


-c.K  J  +  Me.  =  0 
4  4 


From  39,  if  J  r  -  M  *  0,  (i.e.,  K  *  X,  cf.  equation  33) 


c4  ~  0 


From  40  and  36, 


Me  j  -  JK^  cL  =  F 


or  using  34  to  eliminate  M 


Cj  -  F  /  (M  -  JK2)  =  F  /  [J(X2  -  K2)] 


From  37 


c  -  E  /  (M  -  K2J)  =  E  /  [J  (X2-  K2) 


From  38  and  42 


c2  =  [A  (X2  -  K2)  -  2KE]  /  [J  (X2  -  K2)2; 


Therefore  35  simplifies  to 


w  =  c.,sin(Kt  +  9  )  +  c„cos(Kt  +  8  )  +  c_t  sin(Kt  +  9  ) 
vl  0  2  0  3  0 


and  the  complete  solution  is  the  sum  of  44  and  32 

w  =  c,  cost  Xt  +  }>)  +  c  sin(Kt  +  9  )  +  c„cos(Kt  +  9,) 
y  h  1  0  2  0 

+  c3t  sin(Kt  +  9^) 


which  can  be  rewritten  as 


w  =  c  cos(~t)  +  c  sin('t)  +  c  sin(Kt  +  9  )  +  c  cos(Kt  +  9  ) 
y  5  6  i  0  2  0 

-  —  (Ab) 

+  c  t  sin( Kt  +  9  ) 

The  constants  of  integration  tv  and  can  be  evaluated  in  terms  of  the  initial 
conditions  w  and  w  ^  .  Recall  that  c.  ,  c0 ,  and  c,  are  already  determined;  see 
41,  42,  and  IV.  y0  12! 

wy0  =  wy(J[t  =  tc  -  L  /  ( 2R)  ]  =  wy0  [  t  =  0  j  (Ah) 

V  =  >lt  =  V  L/  (2R)  1  =  "y0(t  =  01  (47) 

where  t  is  the  center  time,  the  time  of  the  middle  of  the  burn.  For  conven- 

c 

ience,  make  the  definition 


t  =  L/R  (48) 

which  is  the  burn  time  of  a  thruster.  Thus  the  initial  conditions  become 


£ 

o 

n 

W  (t 
yO 

=  t  -  t/2) 
c 

'  >(t ' 

0) 

(49) 

it 

?. 

•  3 

4-1 

•  3 

=  tc-  t/2) 

*  v1  ■ 

0) 

(50) 

The  derivative 

of  45 

is 

w  ( t) 

y 

=  "XC5 

sin(At)  + 

Ac  cos( At) 
o 

+  Kc^cos(Kt  +  9^) 

-  Kc^  sin( Kt  +  9^)  +  c^sin(Kt  +  0^)  +  c^Kt  cos(Kt  +  9  ) 


(51) 


Thus,  at  t  =  t  -  t/2,  or  t  =  0 
c 


w  =  c  +  c  s i n [  9  ]  +  c  cos (9  ] 
yO  i  1  0  2  0 


(52) 


and 


w 

yo 


Ac  +  Kc  cos[9  ]  -  Kc„sin[sin9  )  +  c  sin[9  ] 
nil)  2  U  i  U 


(53) 


■  -  .  y  -  -  • 
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Solving  for  from  53 


c,  =  1 1 / X ]  {w  -  Kc  cos[ 9  ]  +  Kc  sin[ 9  ]  -  c.sin{0  ]  } 
6  yolo2o3o 


Solving  52  for 


C5  ”  {wy0  '  cisin[0o3  "  c2COsleQ)  } 


Similarly,  the  solution  to  26  is  required 


Jw  +  M  w  =  (Et  +  F)  cos(Kt  +  9,,)  +  A  sin(Kt  +  6  ) 
z  z  0  0 


where 

J  =  I  =*  J  ( I  /I  )  (57) 

z  z  y 

M  =  w2  (I  -  I  )  (I  -  I  )  /  I  =  M  (I  /I  )  (58) 

xxzxy  y  z  y 

E  =>  A  [K  +  w  (I  -  I  )  /  I  ]  (59) 

xx  y  y 

F  =  B  [K  +  w  (I  -  I  )  /I  ]  (60) 

xx  y  y 

The  homogeneous  part  of  56  is  formally  identical  to  the  homogeneous  part  of  27; 
therefore  the  homogeneous  solution  of  56  is  formally  the  same  as  the  homogeneous 
solution  of  27  which  is  expressed  by  32  and  34.  For  the  inhomogeneous  solution, 
again  try  the  form  of  35 

w  =  c  sin( Kt  +  9  )  +  c„cos(Kt  +  9.1  +  c.t  sin(Kt  +  9  ) 
z  1  0  2  0  3  0 

+  c.t  cos(Kt  +  0  )  (81) 

4  0 

Similarly,  differentiate  twice  and  substitute  into  the  differential  equation,  in 
this  case  56,  yielding 

Cj  =  [A(M  -  K  J)  +  2JKF:]  /  [M  -  K2J]2 
2 

or  since  M/J  =  M/J  =  X  ,  then 


m 


c  =  I A( X  -  K  )  +  2KE ]  /  ( J( X  -  K  )  I 


(62) 


c2  =  F/[M  -  K2J]  =  F/lJ(X2  -  K2)] 


=  0 


and 


c4  =  E/  (M  -  K2J)  =  E /  [J(X2  -  K2)  ] 


which  should  be  compared  to  40  through  43.  Following  the  pattern  used 
ate  the  constants  >.i  integration  in  equation  45,  first  write  the  recase 
solution  of  56,  which  is  analogous  to  45,  as 

w  =  crcos(Xt)  +  c.sin(Xt)  +  c, sin(Kt  +  9  )  +  c  cos(Kt  +  9  ) 
z  5  6  1  0  2  0 

+  c^t  cos  (Kt  +  0^) 

where 

C5  =  ^Wz0  "  ^sin^]  -  c2cos[eo)  } 

and 


=  1 1/A]  fw^  -  KCj  cos[ 0^ ]  +  Kc2 sin[ 0^  J  -  c4cos[9QJ  } 

To  obtain  the  final  result,  evaluate  45  and  66  at  t  =  tc  -  t/2  or, 
equivalently,  at  t  =  x  =  L/R. 

w  =  c. cos [ XL/ R]  +  c,  sin(XL/R]  +  c, sin[KL/R  +  9.] 
y  j  b  I  0 

+  c2cos[KL/R  +  Oq]  +  c^ [L/R]  sin[KL/R  +  9^] 

w  =  c  cos(XL/L]  +  cr  sin[XL/R]  +  c ,sin[KL/R  +  9  I 
z  5  6  1  0 

+  c  cos[KL/R  +  0)  +  c [L/R]  cos[KL/R  +  0  ] 


(63) 

(64) 

(65) 

to  evalu- 
complete 

(66) 

(67) 

(68) 


(69) 

(70) 


Care  must  be  taken  in  specifying  initial  conditions.  The  original  equations 
to  be  solved  were  the  coupled  first  order  equations  2  and  3.  For  these  equa¬ 
tions,  only  w  and  w  need  to  be  specified.  However,  when  the  equations  are 

decoupled,  thd  order  is  increased  and  the  initial  conditions,  w  and  w  ,  are 
also  required.  These  cannot  be  specified  arbitrarily.  Instead,  to  avoid  spuri¬ 
ous  solutions  and  to  assure  consistency  with  the  original  coupled  first  order 
equations  2  and  3,  w  and  w  in  principle  must  be  calculated  from  w  and  w 
using  equations  2  a&2  3.  2%quivalently ,  equations  21  and  22  can  ybe  solved 

for  w  and  w  with  t  set  to  zero  to  yield  the  required  initial  conditions,  viz. 
y  z 


w  =  {w  w  [I 
yo  zo  x  z 


w  =  {w  w  [I 
zo  yo  x  x 


I  ]  +  B  cos  9  }  /  I 
x  .  o  y 

I  ]  +  B  sin  9  }  /  I 
y  o  z 


(71) 

(72) 


where  B  could  be  replaced  by  R  IAx.  (See  equation  19.)  A  program  implementing 
this  approach  is  in  the  appendix. 


NUMERICAL  EXAMPLE 


The  analytical  approach  described  above  was  tested  by  comparison  to  numeri¬ 
cal  integration  of  equations  2  and  3  using  a  fourth  order  Runge-Kutta  integration 
technique.  These  results  are  reported  in  table  1.  The  first  entry  is  the  exact 
result  obtained  by  the  analytical  technique.  The  second  entry  represents  numeri¬ 
cal  integration  using  an  integration  time  step  of  0.00000001  second.  Results  are 
interpolated  to  yield  a  value  at  the  exact  time  of  the  end  of  burn.  (If  the 
simulations  are  adjusted  to  make  them  stop  at  the  same  time,  the  numerical  and 
analytical  results  for  w  and  wz  agree  to  at  least  eight  significant  figures.) 
The  third  and  fourth  table  entries  use  integration  time  steps  of  0.0000001  second 
and  0.000001  second,  respectively.  Note  that  the  results  obtained  numerically 
required  from  23  to  2085  times  the  computer  CPU  time  as  the  new  analytical 
approach  to  obtain  3  to  6  significant  figure  accuracy,  respectively. 


CONCLUSIONS 

The  closed  form,  analytical  solution  has  been  accurate  and  highly  efficient 
for  the  simulation  of  the  flight  of  projectiles  with  maneuvers  implemented  by 
helical,  quas i-impul si ve  thrusters.  The  technique  will  be  implemented  in  the 
further  development  of  such  projectiles  at  ARDEC. 


Table  1.  Comparison  of  results 


Model  type 

Time  step 
(sec) 

CPU  time 
(sec) 

v 

(rad) 

(rad) 

Analytical 

N/A 

0.0106 

0.53275603 

-0.41630908 

Numerical 

10-8 

22.096 

0.53275601 

-0.41630910 

Numerical 

10"7 

2.278 

0.53275061 

-0.41631581 

Numerical 

10-6 

0.246 

0.53234719 

-0.41681537 

SYMBOLS 

'"‘-1 

F 

Force  applied  to  body  by  burning  propellant  (pounds) 

;v; 

D 

Diameter  of  cylindrical  body  section  (feet) 

$ 

I 

Specific  impulse  of  thruster  (pound-seconds/foot  of  thruster) 

:»s 

*x’*y’*z 

Moments  of  inertia  about  fixed  axes  of  body  (slug-feet  squared) 

i’j 

vfl 

1 

Distance  down  groove  (feet) 

ii 

•'ll* 

L 

Total  length  of  thruster  groove  (feet) 

N  N  N 
x  ’  Ny  ,INz 

Components  of  torque  applied  to  body  (pound-feet) 

*lu 

it 

R 

Linear  burning  rate  of  propellant  (feet/second) 

S 

t 

Time  (seconds) 

•$£ 

t 

Time  measured  from  zero  at  start  of  burn  (seconds) 

ft 

tc 

Time  t  at  center  of  thruster  burn 

$ 

T 

Twist  of  helix  (turns  per  caliber) 

*3 

x,y 

Coordinates  of  point  on  groove  when  "unwrapped"  (feet) 

$ 

x,y,z 

Coordinates  of  point  in  body  frame  (feet) 

1 

0 

Polar  angle  of  point  on  groove,  positive  from  Y  toward  Z  (radians) 

yi 

T 

Time  of  a  thruster  burn  (seconds) 

$L 

$ 

\ 

Natural  frequency  of  solution  to  homogeneous  differential  equation 
(per  second) 

wx’VwZ 

Components  of  spin  in  body  frame  (radians  per  second) 

;,yi 

& 

Mi 

$ 

vv 

i*7 

•;*; 

;> 

•  ■» 

• 

*\t 
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PROGRAM  CGS P 


17 


onnnnonn  non  o  n  r. 


THIS  IS  A  SAMPLE  DRIVER  PROGRAM  FOR  TESTING  SUBROUTINE  OMEGA 

COMMON  R ,  L , I ,  T ,  D , WX , IX,IY,IZ,PI , WYZ ERO , WZZ ERO , WYZDOT , WZZDOT , 
+  THZER, XZERO,FREQN,FREQF 
REAL  L, I , IX, IY, IZ 

DEFINITIONS 

R=2 . 23! 0E04 
L  =  1 . 3  3  3 
1=12. 842 
T=-l . 5755E-02 
D= . 4 199 
WX=1000 . 0 
IX=. 136 
IY=. 930 
I  Z  =  I  Y 

PI=3. 141392654 

BURN  OUT  TIME.  THE  FOLLOWING  EXACT  VALUE  MAY  BE  REPLACED 
BY  THE  ACTUAL  INTEGRATION  CUT  OFF  TIME  IN  A  NUMERICAL 
INTEGRATION  ROUTINE  FOR  PURPOSES  OF  COMPARISON. 

TIME  =  L/R 

INITIAL  CONDITIONS 


C 


C 


C 

c 

c 


WYZERO  =0.1 
WZZERO  =  0.2 

XZ  ERO  =  0.6667 
THZER  =  PI/7. 

CALL  OMEGA (TIME, WY,WZ) 

PRINT*,'  TIME  WY  WZ ' 

WRITE  (  * , 100 )  TIME, WY , WZ 
WRITE  (* ,101) 

1  CONTINUE 
STOP 

100  FORMAT (4F15. 12,//) 

101  FORMAT ( ' 1 ' ) 

END 

SUBROUTINE  OMEGA (T IME , WY , WZ ) 

COMMON  R,L,I,T,D, WX ,IX,IY,IZ,PI , WYZ ERO, WZZERO, WYZDOT, WZZDOT 
+  , THZER, XZERO,FREQN,FREQF 

REAL  L, I, IX, IY, IZ, J,JBAR 

FREQN  =  NATHP^L  FREQUENCY  OF  THE  SYSTEM 
FREQF  =  FORCING  FUNCTION  FREQUENCY 
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non  noon  non 


EL=SQRT(1.  +  (PI *T ) * ( PI *T ) ) 

FREQN=SQRT (  (WX*WX ) * ( I X  — I Z ) * (IX-IY) / ( I Y* I Z )  ) 

FREQF  =  -(2.  *  PI  *  T  *  R  )/(  D*EL  ) 

J  =  I Y 
JBAR= I Z 

A  =  -(R*R)  *  I/EL 

B=  R  *  I  *  XZERO 

E=  -A  *  ((  W X  *  (IX-IZ)/IZ)  +  FREQF) 

EBAR=  + A  *  ( (  WX  *  ( IX-IY ) /I Y )  +  FREQF) 

F=  -B  *  ( (  WX  *  (IX-IZ)/IZ)  +  FREQF) 

FBAR=  B  *  {  (  WX  *  ( I X-I Y ) /I Y )  +  FRFQ^) 

DENOM=  (FREQN*FREQN )  -  (FREQF*FREQF) 

CONSTANTS  DERIVED  FROM  THE  DECOUPLING  PROCESS. 

C1=F/ ( J  *DENOM) 

C2B=FBAR/ ( JBAR*DENOM) 

C2=  (A/(J*DENOM)  )-  (E*2 . * FREQF/ ( J* (DENOM*DENOM)  )  ) 

C1B= (A/(JBAR*DENOM)  ) + (EBAR*2 . *FREQF/ ( JBAR* (DENOM*DENOM)  )  ) 

C3  =  E/ ( J  *DENOM) 

C4B=EBAR/ ( JBAR*DENOM) 

C4  =  0 . 

C3B=0 . 

SUPPRESS  SPURIOUS  SOLUTIONS  INTRODUCED  BY 
TAKING  SECOND  DERIVATIVES  OF  WY  AND  WZ 

WYZDOT  =  (WZZERO*WX* (IZ-IX)  +  B*COS (THZER) )/IY 
WZZDOT  =  (WYZERO*WX*  (IX-IY)  +  B*S I N  (THZER )  ) /I Z 

CONSTANTS  DERIVED  FROM  THE  INITIAL  CONDITIONS 

C5=  (WYZERO-Cl*SIN (THZER) -C2*COS  (THZER)  ) 

C6  =  (1. /FREQN) * (WYZDOT-FREQF*Cl*COS  (THZER) 

+  +  FREQF*C2*SIN  (THZER)  -  C3*S IN (THZER)  ) 

C5B=  (WZZER0-C1B*SIN (THZER)  -C2B*COS (THZER)  ) 

C6B  =  (1 ,/FREQN) * (WZZDOT 

+  -FREQF*C1B*C0S  (THZER) +FREQF* C2B* S I N (THZER) -C4B*COS (THZER)  ) 
C 

WY  =  C5*COS  (FREQN*TIME)  +  C6 * S I N ( FREQN *T IME ) 

+  +  C1*SIN(FREQF*TIME+THZER)  + 

+  C2*COS ( FREQF*T IME  +  THZ ER )  +  C 3 *T I ME* S I N  ( FREQF*T I ME  +  THZ ER ) 

C 

WZ  =  C5B*COS (FREQN*TIME)  +  C6B*  S I N ( FREQN  *T IME) 

+  +  C 1B*S I N (FREQF*T IME  +  THZ  ER )  +  C 2 B*COS ( FREQF*T I ME  +  THZER ) 

+  +C4B*TIME*COS ( FREQF* T I ME+THZER) 

C 

RETURN 

END 


Sample  Execution 


TIME 

0.0000;"  9748991 


WY  WZ 

0.532756032730  -0.416309075410 
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